#################################################
####       Replication Code for              ####
####  A Re-Assessment of Reporting Bias in   ####
#### Event-Based Violence Data with Respect  ####
####         to Cell Phone Coverage          ####
#################################################
####                                         ####  
####                1/18/2017                ####         
####                 Generate all Figures    ####  
#################################################


library(haven)
library(foreign)
library(ggplot2)
library(plyr)
library(coda)

# set working directory to the replication folder


##########################################
##### Main Paper
##########################################


## Create Figure 1: 

# panel (a)
load("results_Windows_high.rda")
mean = NULL
lower = NULL
upper = NULL
window = NULL
for(i in 1:45){
  mean = c(mean, results_Windows_high[[i]][1])
  lower = c(lower, results_Windows_high[[i]][2])
  upper = c(upper, results_Windows_high[[i]][3])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_high.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Windows", y="Coefficient Estimate") 
plot(coef_plot)
dev.off()



# b)
load("results_SpatWindows.rda")

# panel (b) 
# this loads the results from the window analysis
# this part is to estimate the comparison risk ratio from the full data
#load("EventsPerWindow.rda")
#load("africa_shape_grid_windows.RData")
#load("W_windows.Rdata")

mean = NULL
lower = NULL
upper = NULL
window = NULL
for(i in 1:45){
  mean = c(mean, results_SpatWindows[[i]][1])
  lower = c(lower, results_SpatWindows[[i]][2])
  upper = c(upper, results_SpatWindows[[i]][3])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_spatialWindows.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Windows", y="Coefficient Estimate") 
plot(coef_plot)
dev.off()



# panel (c) 
win <- read_dta("windows7_1.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_high.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimate")
plot(coef_plot)
dev.off()



# panel (d) 
win <- read_dta("windows1.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimate")
plot(coef_plot)
dev.off()

## Create Figure 2: 

# panels (a)-(e)
load("simu_out_CS_predicted.rda")

simu_out_CS_predicted = data.frame(simu_out_CS_predicted)
names(simu_out_CS_predicted) = c("iter","added","coef","SE")
unique(simu_out_CS_predicted$added)

results95 <- simu_out_CS_predicted[simu_out_CS_predicted$added==10,]
results90 <- simu_out_CS_predicted[simu_out_CS_predicted$added==20,]
results85 <- simu_out_CS_predicted[simu_out_CS_predicted$added==32,]
results80 <- simu_out_CS_predicted[simu_out_CS_predicted$added==46,]
results75 <- simu_out_CS_predicted[simu_out_CS_predicted$added==61,]
results70 <- simu_out_CS_predicted[simu_out_CS_predicted$added==78,]
results65 <- simu_out_CS_predicted[simu_out_CS_predicted$added==98,]


draws95 <- NA 
draws90 <- NA 
draws85 <- NA 
draws80 <- NA 
draws75 <- NA
draws70 <- NA
draws65 <- NA

for (i in 1:1000){
  draws95 <- c(draws95,rnorm(200,mean=results95[i,"coef"],sd=results95[i,"SE"]))
  draws90 <- c(draws90,rnorm(200,mean=results90[i,"coef"],sd=results90[i,"SE"]))
  draws85 <- c(draws85,rnorm(200,mean=results85[i,"coef"],sd=results85[i,"SE"]))
  draws80 <- c(draws80,rnorm(200,mean=results80[i,"coef"],sd=results80[i,"SE"]))
  draws75 <- c(draws75,rnorm(200,mean=results75[i,"coef"],sd=results75[i,"SE"]))
  draws70 <- c(draws70,rnorm(200,mean=results70[i,"coef"],sd=results70[i,"SE"]))
  draws65 <- c(draws65,rnorm(200,mean=results65[i,"coef"],sd=results65[i,"SE"]))
  
}

draws95 <- draws95[-1] 
draws90 <- draws90[-1] 
draws85 <- draws85[-1] 
draws80 <- draws80[-1] 
draws75 <- draws75[-1] 
draws70 <- draws70[-1] 
draws65 <- draws65[-1] 


#plots that include the estimation uncertainty

#95%
greater0 <- ifelse(draws95>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws95)$x, y=density(draws95)$y)
summary(gg)
# generate the plot
pdf("coef_95_CS_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.009, y =75, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#90%
greater0 <- ifelse(draws90>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws90)$x, y=density(draws90)$y)
summary(gg)
# generate the plot
pdf("coef_90_CS_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.009, y =75, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#85%
greater0 <- ifelse(draws85>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws85)$x, y=density(draws85)$y)
summary(gg)
# generate the plot
pdf("coef_85_CS_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.008, y =65, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()


#80%
greater0 <- ifelse(draws80>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws80)$x, y=density(draws80)$y)
summary(gg)
# generate the plot
pdf("coef_80_CS_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.009, y =65, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#75%
greater0 <- ifelse(draws75>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws75)$x, y=density(draws75)$y)
summary(gg)
# generate the plot
pdf("coef_75_CS_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.033, y =60, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()


################################################
#### Appendix
################################################







## Create Figure 1
# panel (a)
# this loads the results from the window analysis
load("results_windows_cmean_high.rda")
mean = NULL
lower = NULL
upper = NULL
window = NULL
for(i in 1:45){
  mean = c(mean, results_windows[[i]][1])
  lower = c(lower, results_windows[[i]][2])
  upper = c(upper, results_windows[[i]][3])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_cmean_high.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates") #+ scale_x_continuous(breaks = c(1,5,10,15,20,25,30,35,40,45), labels = round(EventsPerWindow$win2008[c(1,5,10,15,20,25,30,35,40,45)],2)) 
plot(coef_plot)
dev.off()


# panel (b)
win <- read_dta("windows12.dta")
win <- subset(win,win$parm=="cell_dum_07")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_cross_fe_high_spat.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()




###windowed random and with best/high ordering

# panel (a)
load("random_high_results.rda")
names(random_high_results) = c("window", "iteration", "EventsInWindow", "coefficients")

random_high_results = random_high_results[order(random_high_results$window, random_high_results$iteration),]


data = ddply(random_high_results,.variable = "window",summarise, mean = mean(coefficients), lower = quantile(coefficients, 0.025), upper = quantile(coefficients, 0.9725))
pdf("Coef_Windows_random_high.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Windows", y="Coefficient Estimate") 
plot(coef_plot)
dev.off()

##panel (b) 
load("random_best_results.rda")

names(random_best_results) = c("window", "iteration", "EventsInWindow", "coefficients")
random_best_results = random_best_results[order(random_best_results$window, random_best_results$iteration),]


data = ddply(random_best_results,.variable = "window",summarise, mean = mean(coefficients), lower = quantile(coefficients, 0.025), upper = quantile(coefficients, 0.9725))
pdf("Coef_Windows_random_best.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Windows", y="Coefficient Estimate") 
plot(coef_plot)
dev.off()



# figure 3
load("results_SpatWindows_BestHigh.rda")
mean = NULL
lower = NULL
upper = NULL
window = NULL
for(i in 1:45){
  mean = c(mean, results_SpatWindows_BestHigh[[i]][1])
  lower = c(lower, results_SpatWindows_BestHigh[[i]][2])
  upper = c(upper, results_SpatWindows_BestHigh[[i]][3])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_besthigh.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) 
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates") #+ scale_x_continuous(breaks = c(1,5,10,15,20,25,30,35,40,45), labels = round(EventsPerWindow$win2008[c(1,5,10,15,20,25,30,35,40,45)],2)) 
plot(coef_plot)
dev.off()






## Create Figure 2
# panel (a)
win <- read_dta("windows2.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_prec.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()


# panel (b)
win <- read_dta("windows5.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_new.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()



## Create Figure 3
# panel (a)
win <- read_dta("windows7.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_high_prec.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()


# panel (b)
win <- read_dta("windows10.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_high_new.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()


## Create Figure 4
win <- read_dta("windows3.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_count.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()

## Create Figure 5
win <- read_dta("windows8.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_high_count.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()

## Figure 6
# panel (a)
win <- read_dta("windows4.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_count_prec.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()

#panel (b)
win <- read_dta("windows6.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_count_new.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()


## Figure 7
# panel (a)
win <- read_dta("windows9.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_high_prec_count.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()

#panel (b)
win <- read_dta("windows11.dta")
win <- subset(win,win$parm=="cell")

mean = NULL
lower = NULL
upper = NULL
window = NULL

for(i in 1:45){
  mean = c(mean, win$estimate[i])
  lower = c(lower, win$min95[i])
  upper = c(upper, win$max95[i])
  window = c(window,i)
}
data = data.frame(window,mean,lower,upper)

pdf("Coef_Windows_panel_high_new_count.pdf",height=7,width=7)
coef_plot = ggplot(data)
coef_plot = coef_plot + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
coef_plot = coef_plot + geom_linerange(aes(x = window, ymin = lower,ymax = upper),lwd = 1)
coef_plot = coef_plot + geom_point(aes(x = window, y = mean),lwd = 3,shape = 21, fill = "WHITE")
coef_plot = coef_plot  + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=15)) + theme(axis.text.y=element_text(size=15))
coef_plot <- coef_plot + theme(axis.title.x = element_text(size = 20),axis.title.y = element_text(size = 20)) + labs(title = "", x="Windows", y="Coefficient Estimates")
plot(coef_plot)
dev.off()



## Figure 8
load("simu_out_Panel_predicted.rda")
simu_out_Panel_predicted = data.frame(simu_out_Panel_predicted)
names(simu_out_Panel_predicted) = c("iteration","events_added","coef","SE")

results95 <- simu_out_Panel_predicted[simu_out_Panel_predicted$events_added==28,]
results90 <- simu_out_Panel_predicted[simu_out_Panel_predicted$events_added==59,]
results85 <- simu_out_Panel_predicted[simu_out_Panel_predicted$events_added==93,]
results80 <- simu_out_Panel_predicted[simu_out_Panel_predicted$events_added==132,]
results75 <- simu_out_Panel_predicted[simu_out_Panel_predicted$events_added==176,]

draws95 <- NA 
draws90 <- NA 
draws85 <- NA 
draws80 <- NA 
draws75 <- NA 

for (i in 1:1000){
  draws95 <- c(draws95,rnorm(200,mean=results95[i,"coef"],sd=results95[i,"SE"]))
  draws90 <- c(draws90,rnorm(200,mean=results90[i,"coef"],sd=results90[i,"SE"]))
  draws85 <- c(draws85,rnorm(200,mean=results85[i,"coef"],sd=results85[i,"SE"]))
  draws80 <- c(draws80,rnorm(200,mean=results80[i,"coef"],sd=results80[i,"SE"]))
  draws75 <- c(draws75,rnorm(200,mean=results75[i,"coef"],sd=results75[i,"SE"]))
}

draws95 <- draws95[-1] 
draws90 <- draws90[-1] 
draws85 <- draws85[-1] 
draws80 <- draws80[-1] 
draws75 <- draws75[-1] 


#plots that include the estimation uncertainty

#95%
greater0 <- ifelse(draws95>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws95)$x, y=density(draws95)$y)
summary(gg)
# generate the plot
pdf("coef_95_panel_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.025, y =75, label = "99.8% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#90%
greater0 <- ifelse(draws90>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws90)$x, y=density(draws90)$y)
summary(gg)
# generate the plot
pdf("coef_90_panel_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.025, y =75, label = "99% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#85%
greater0 <- ifelse(draws85>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws85)$x, y=density(draws85)$y)
summary(gg)
# generate the plot
pdf("coef_85_panel_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.02, y =65, label = "97% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()


#80%
greater0 <- ifelse(draws80>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws80)$x, y=density(draws80)$y)
summary(gg)
# generate the plot
pdf("coef_80_panel_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.02, y =65, label = "93% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#75%
greater0 <- ifelse(draws75>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws75)$x, y=density(draws75)$y)
summary(gg)
# generate the plot
pdf("coef_75_panel_predicted.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.02, y =60, label = "86% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()




## Figure 9

load("simu_out_CS_completeRandom.rda")
simu_out_CS_completeRandom = data.frame(simu_out_CS_completeRandom)
names(simu_out_CS_completeRandom) = c("prop","coef","SE")

results95 <- simu_out_CS_completeRandom[simu_out_CS_completeRandom$prop==0.95,]
results90 <- simu_out_CS_completeRandom[simu_out_CS_completeRandom$prop==0.9,]
results85 <- simu_out_CS_completeRandom[simu_out_CS_completeRandom$prop==0.85,]
results80 <- simu_out_CS_completeRandom[simu_out_CS_completeRandom$prop==0.8,]
results75 <- simu_out_CS_completeRandom[simu_out_CS_completeRandom$prop==0.75,]
results70 <- simu_out_CS_completeRandom[simu_out_CS_completeRandom$prop==0.7,]


draws95 <- NA 
draws90 <- NA 
draws85 <- NA 
draws80 <- NA 
draws75 <- NA
draws70 <- NA

for (i in 1:1000){
  draws95 <- c(draws95,rnorm(200,mean=results95[i,"coef"],sd=results95[i,"SE"]))
  draws90 <- c(draws90,rnorm(200,mean=results90[i,"coef"],sd=results90[i,"SE"]))
  draws85 <- c(draws85,rnorm(200,mean=results85[i,"coef"],sd=results85[i,"SE"]))
  draws80 <- c(draws80,rnorm(200,mean=results80[i,"coef"],sd=results80[i,"SE"]))
  draws75 <- c(draws75,rnorm(200,mean=results75[i,"coef"],sd=results75[i,"SE"]))
  draws70 <- c(draws70,rnorm(200,mean=results70[i,"coef"],sd=results70[i,"SE"]))
  
}

draws95 <- draws95[-1] 
draws90 <- draws90[-1] 
draws85 <- draws85[-1] 
draws80 <- draws80[-1] 
draws75 <- draws75[-1] 
draws70 <- draws70[-1] 


#plots that include the estimation uncertainty

#95%
greater0 <- ifelse(draws95>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws95)$x, y=density(draws95)$y)
summary(gg)
# generate the plot
pdf("coef_95_CS_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.009, y =75, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#90%
greater0 <- ifelse(draws90>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws90)$x, y=density(draws90)$y)
summary(gg)
# generate the plot
pdf("coef_90_CS_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.009, y =75, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#85%
greater0 <- ifelse(draws85>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws85)$x, y=density(draws85)$y)
summary(gg)
# generate the plot
pdf("coef_85_CS_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.008, y =65, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()


#80%
greater0 <- ifelse(draws80>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws80)$x, y=density(draws80)$y)
summary(gg)
# generate the plot
pdf("coef_80_CS_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.009, y =65, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#75%
greater0 <- ifelse(draws75>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws75)$x, y=density(draws75)$y)
summary(gg)
# generate the plot
pdf("coef_75_CS_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred")
p = p + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.033, y =60, label = "100% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

## Figure 10
load("simu_out_Panel_completeRandom.rda")

names(simu_out_Panel_completeRandom) = c("iteration","events_added","coef","SE","pval")

results95 <- simu_out_Panel_completeRandom[simu_out_Panel_completeRandom$events_added==28,]
results90 <- simu_out_Panel_completeRandom[simu_out_Panel_completeRandom$events_added==59,]
results85 <- simu_out_Panel_completeRandom[simu_out_Panel_completeRandom$events_added==93,]
results80 <- simu_out_Panel_completeRandom[simu_out_Panel_completeRandom$events_added==132,]
results75 <- simu_out_Panel_completeRandom[simu_out_Panel_completeRandom$events_added==176,]

draws95 <- NA 
draws90 <- NA 
draws85 <- NA 
draws80 <- NA 
draws75 <- NA 

for (i in 1:1000){
  draws95 <- c(draws95,rnorm(200,mean=results95[i,"coef"],sd=results95[i,"SE"]))
  draws90 <- c(draws90,rnorm(200,mean=results90[i,"coef"],sd=results90[i,"SE"]))
  draws85 <- c(draws85,rnorm(200,mean=results85[i,"coef"],sd=results85[i,"SE"]))
  draws80 <- c(draws80,rnorm(200,mean=results80[i,"coef"],sd=results80[i,"SE"]))
  draws75 <- c(draws75,rnorm(200,mean=results75[i,"coef"],sd=results75[i,"SE"]))
}

draws95 <- draws95[-1] 
draws90 <- draws90[-1] 
draws85 <- draws85[-1] 
draws80 <- draws80[-1] 
draws75 <- draws75[-1] 


#plots that include the estimation uncertainty

#95%
greater0 <- ifelse(draws95>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws95)$x, y=density(draws95)$y)
summary(gg)
# generate the plot
pdf("coef_95_Panel_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.025, y =75, label = "98% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#90%
greater0 <- ifelse(draws90>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws90)$x, y=density(draws90)$y)
summary(gg)
# generate the plot
pdf("coef_90_Panel_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.025, y =75, label = "95% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#85%
greater0 <- ifelse(draws85>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws85)$x, y=density(draws85)$y)
summary(gg)
# generate the plot
pdf("coef_85_Panel_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.02, y =65, label = "89% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()


#80%
greater0 <- ifelse(draws80>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws80)$x, y=density(draws80)$y)
summary(gg)
# generate the plot
pdf("coef_80_Panel_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.02, y =65, label = "79% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()

#75%
greater0 <- ifelse(draws75>0,1,0)
mean(greater0)
#  generate kdf
gg <- data.frame(x=density(draws75)$x, y=density(draws75)$y)
summary(gg)
# generate the plot
pdf("coef_75_Panel_completeRandom.pdf",height=7,width=7)
p =  ggplot(gg,aes(x=x)) + geom_line(aes(x=x,y=y))  + geom_ribbon(data=subset(gg,x>0),aes(x=x,ymin=0,ymax=y))  + geom_vline(xintercept = 0, colour="darkred") + theme(panel.background = element_rect(fill='white', colour='gray')) + theme(legend.position="none") + theme(axis.text.x=element_text(size=12)) + theme(axis.text.y=element_text(size=12)) + annotate("text", x = 0.02, y =60, label = "66% of coefficients \n greater than zero")
p = p + theme(axis.title.x = element_text(size = 15),axis.title.y = element_text(size = 15)) + labs(title = "", x="Coefficient Estimates", y="Density") 
plot(p)
dev.off()




















